GARCH-EVT-COPULA APPROACH (FINAL 2)

IMPORT ALL LIBRARIES

#Import all libraries
library('quantmod')

## Loading required package: xts

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## Loading required package: TTR

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

## Version 0.4-0 included new data defaults. See ?getSymbols.

library('tidyquant')

## Warning: package 'tidyquant' was built under R version 4.0.3

## Loading required package: lubridate

## Warning: package 'lubridate' was built under R version 4.0.3

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

## Loading required package: PerformanceAnalytics

## Warning: package 'PerformanceAnalytics' was built under R version 4.0.3

## 
## Attaching package: 'PerformanceAnalytics'

## The following object is masked from 'package:graphics':
## 
##     legend

## == Need to Learn tidyquant? =====================================================================
## Business Science offers a 1-hour course - Learning Lab #9: Performance Analysis & Portfolio Optimization with tidyquant!
## </> Learn more at: https://university.business-science.io/p/learning-labs-pro </>

library('rugarch')

## Warning: package 'rugarch' was built under R version 4.0.3

## Loading required package: parallel

## 
## Attaching package: 'rugarch'

## The following object is masked from 'package:stats':
## 
##     sigma

library('copula')

## Warning: package 'copula' was built under R version 4.0.3

## 
## Attaching package: 'copula'

## The following object is masked from 'package:lubridate':
## 
##     interval

library('VGAM')

## Warning: package 'VGAM' was built under R version 4.0.3

## Loading required package: stats4

## Loading required package: splines

## 
## Attaching package: 'VGAM'

## The following objects are masked from 'package:copula':
## 
##     log1mexp, log1pexp, rlog

library('ggplot2')
library('GoFKernel')

## Warning: package 'GoFKernel' was built under R version 4.0.3

## Loading required package: KernSmooth

## KernSmooth 2.23 loaded
## Copyright M. P. Wand 1997-2009

library('mistr')

## Warning: package 'mistr' was built under R version 4.0.3

## 
## Attaching package: 'mistr'

## The following objects are masked from 'package:VGAM':
## 
##     dfrechet, dgumbel, dpareto, pfrechet, pgumbel, ppareto, qfrechet,
##     qgumbel, qpareto, rfrechet, rgumbel, rpareto

## The following object is masked from 'package:base':
## 
##     q

library('tseries')

## Warning: package 'tseries' was built under R version 4.0.3

library('stats')
library('fDMA')

## Warning: package 'fDMA' was built under R version 4.0.3

INPUT AND CLEAN THE DATA

CTG<-read.csv('C:/Users/Win 10/Documents/excel_ctg.csv')
MSN<-read.csv('C:/Users/Win 10/Documents/excel_msn.csv')
VIC<-read.csv('C:/Users/Win 10/Documents/excel_vic.csv')
VNM<-read.csv('C:/Users/Win 10/Documents/excel_vnm.csv')


cleandata<-function(df)
  {
  names(df)<-c('Ticker','Date','Open','High','Low','Close','Volume')
  df<-df[,c(2,6)]
  df <- transform(df, Date = as.Date(as.character(Date), "%Y%m%d"))
  df<-df[order(as.Date(df$Date, format="%Y%m%d")),]
  df<-xts(x =df$Close, order.by =df$Date)
  df
}

CTG<-cleandata(CTG)
MSN<-cleandata(MSN)
VIC<-cleandata(VIC)
VNM<-cleandata(VNM)

# Length of each stock is 1733
CTG<-CTG["2013-01-01/2019-12-16"]
MSN<-MSN["2013-01-01/2019-12-16"]
VIC<-VIC["2013-01-01/2019-12-16"]
VNM<-VNM["2013-01-01/2019-12-16"]


# plot the stock price
layout(matrix(c(1,2,3,4),2, 2, byrow = TRUE))
plot(CTG,col="darkorange2",lwd=1.5)
plot(MSN,col="deeppink3",lwd=1.5)
plot(VIC,col="red",lwd=1.5)
plot(VNM,col="blue",lwd=1.5)

layout(matrix(c(1),1,1))
plot(cbind(CTG,MSN,VIC,VNM),col=c("darkorange2","deeppink3","red","blue"),lwd=1.5)

GET THE LOG RETURN

#Calculate the log differences price 
log_diff<-function(x){
    data=log(x)
    data_diff=diff(data,differences = 1)
    data_diff=data_diff[-1,]
    return(data_diff)
}

#Get the log-return for each company
CTG_log_diff<-log_diff(CTG)
MSN_log_diff<-log_diff(MSN)
VIC_log_diff<-log_diff(VIC)
VNM_log_diff<-log_diff(VNM)


#Plot the histogram
layout(matrix(c(1,2,3,4), 2, 2, byrow = TRUE))
hist(CTG_log_diff,breaks=50,col="darkorange2")
hist(MSN_log_diff,breaks=50,col="deeppink3")
hist(VIC_log_diff,breaks=50,col="red")
hist(VNM_log_diff,breaks=50,col="blue")

#Plot the qq-norm and pp-line
layout(matrix(c(1,2,3,4), 2, 2, byrow = TRUE))

qqnorm(CTG_log_diff,main='CTG log diff')
qqline(CTG_log_diff,col="darkorange2")

qqnorm(MSN_log_diff,main='MSN log diff')
qqline(MSN_log_diff,col="deeppink3")

qqnorm(VIC_log_diff,main='VIC log diff')
qqline(VIC_log_diff,col="red")

qqnorm(VNM_log_diff,main='VNM log diff')
qqline(VNM_log_diff,col="blue")

#Summary of the log-return for each company
CTG_sumary<-c(mean(CTG_log_diff),sd(CTG_log_diff),skewness(CTG_log_diff),kurtosis(CTG_log_diff))

MSN_sumary<-c(mean(MSN_log_diff),sd(MSN_log_diff),skewness(MSN_log_diff),kurtosis(MSN_log_diff))

VIC_sumary<-c(mean(VIC_log_diff),sd(VIC_log_diff),skewness(VIC_log_diff),kurtosis(VIC_log_diff))

VNM_sumary<-c(mean(VNM_log_diff),sd(VNM_log_diff),skewness(VNM_log_diff),kurtosis(VNM_log_diff))


sumary_return<-cbind(CTG_sumary,MSN_sumary,VIC_sumary,VNM_sumary)
dimnames(sumary_return)<-list(c("mean","sd","skewness","kurtosis"),c("CTG","MSN","VIC","VNM"))
sumary_return

##                   CTG           MSN         VIC         VNM
## mean     0.0001712405 -7.208537e-05 0.001051504 0.000731781
## sd       0.0192412730  1.884509e-02 0.016515991 0.014155363
## skewness 0.1823632848  9.649796e-02 0.188060157 0.249918131
## kurtosis 2.8227976619  2.316779e+00 3.469084790 2.393312625

#Plot the log returns of each company
layout(matrix(c(1,2,3,4),2,2,byrow=TRUE))
plot(CTG_log_diff,main='CTG log returns',col="darkorange2",lwd=1.5)
plot(MSN_log_diff,main='MSN log returns',col="deeppink3",lwd=1.5)
plot(VIC_log_diff,main='VIC log returns',col="red",lwd=1.5)
plot(VNM_log_diff,main='VNM log returns',col="blue",lwd=1.5)

#Plot the ACF and PACF (auto-correlation function)
layout(matrix(c(1,2,3,4), 2, 2, byrow = TRUE))

acf(CTG_log_diff,main="ACF (RETURN-CTG)")
acf(MSN_log_diff,main="ACF (RETURN-MSN)")
acf(VIC_log_diff,main="ACF (RETURN-VIC)")
acf(VNM_log_diff,main="ACF (RETURN-VNM)")

pacf(CTG_log_diff,main="PACF (RETURN-CTG)")
pacf(MSN_log_diff,main="PACF (RETURN-CTG)")
pacf(VIC_log_diff,main="PACF (RETURN-CTG)")
pacf(VNM_log_diff,main="PACF (RETURN-CTG)")

FIND THE BEST ARMA ORDER WITH SMALLEST AIC (p,q= 0,1,2,3,4)

# Create a function to find order p,q for ARMA (p,q) with corresponding AIC
## This function will return order of ARIMA and its AIC in 1 vector
order_ARMA<-function(p_max,q_max,data)
{
  final.aic <- Inf
final.order <- c(0,0,0)
for (p in 0:p_max) for (q in 0:q_max) {
  if ( p == 0 && q == 0) {
    next
  }
  
  arimaFit = tryCatch( arima(data, order=c(p, 0, q)),
                       error=function( err ) FALSE,
                       warning=function( err ) FALSE )
  
  if( !is.logical( arimaFit ) ) {
    current.aic <- AIC(arimaFit)
    if (current.aic < final.aic) {
      final.aic <- current.aic
      final.order <- c(p, 0, q)
      final.arima <- arima(data, order=final.order)
    }
  } else {
    next
  }
}
return(c(final.order,final.aic))
}

#Find order for each company
order_CTG<-order_ARMA(4,4,CTG_log_diff)
order_CTG

## [1]     1.000     0.000     1.000 -8788.871

order_MSN<-order_ARMA(2,2,MSN_log_diff)
order_MSN

## [1]     2.000     0.000     2.000 -8848.717

order_VIC<-order_ARMA(4,4,VIC_log_diff)
order_VIC

## [1]     0.00     0.00     1.00 -9294.27

order_VNM<-order_ARMA(3,3,VNM_log_diff)
order_VNM

## [1]     3.000     0.000     3.000 -9833.654

FIT THE ARMA-GARCH MODEL

CTG COMPANY

#Fit the garch model to CTG
garchspec_CTG <- ugarchspec(
  mean.model=list(armaOrder=c(order_CTG[1],order_CTG[3])),
  variance.model=list(model="gjrGARCH",garchOrder = c(1, 1)),
  distribution.model = "sstd")
garchfit_CTG<-ugarchfit(data=CTG_log_diff,spec=garchspec_CTG)
coef(garchfit_CTG)

##            mu           ar1           ma1         omega        alpha1 
## -4.358787e-04 -9.229337e-01  9.119737e-01  7.922143e-06  9.143981e-02 
##         beta1        gamma1          skew         shape 
##  8.497311e-01  1.096070e-01  1.129772e+00  4.471015e+00

#Standardized residuals for CTG
standardize_residual_CTG<-residuals(garchfit_CTG)/sigma(garchfit_CTG)

# Plotting
layout(matrix(c(1,2,3,4),2,2,byrow=TRUE))
acf(standardize_residual_CTG,main="ACF (RESIDUAL-CTG)")
pacf(standardize_residual_CTG, main="PACF (RESIDUAL-CTG)")

plot(standardize_residual_CTG,main="RESIDUAL-CTG",col="darkorange2",lwd=1)
qqnorm(standardize_residual_CTG,main="QQ-PLOT CTG")
qqline(standardize_residual_CTG,col="darkorange2")

MSN COMPANY

#Fit the garch model to MSN
garchspec_MSN <- ugarchspec(
  mean.model=list(armaOrder=c(order_MSN[1],order_MSN[3])),
  variance.model=list(model="gjrGARCH",garchOrder = c(1, 1)),
  distribution.model = "sstd")
garchfit_MSN<-ugarchfit(data=MSN_log_diff,spec=garchspec_MSN)
coef(garchfit_MSN)

##            mu           ar1           ar2           ma1           ma2 
## -4.747141e-04 -2.264270e-02 -8.158039e-01  2.341002e-02  7.747639e-01 
##         omega        alpha1         beta1        gamma1          skew 
##  1.661256e-05  1.674318e-01  8.251220e-01  1.271868e-02  1.023722e+00 
##         shape 
##  3.245284e+00

#Standardized residuals for MSN
standardize_residual_MSN<-residuals(garchfit_MSN)/sigma(garchfit_MSN)

# plotting
layout(matrix(c(1,2,3,4),2,2,byrow=TRUE))
acf(standardize_residual_MSN,main="ACF (RESIDUAL-MSN)")
pacf(standardize_residual_MSN,main="PACF (RESIDUAL-MSN)")

plot(standardize_residual_MSN,col="deeppink3",main="RESIDUAL-MSN",lwd=1)
qqnorm(standardize_residual_MSN,main="QQ-PLOT CTG")
qqline(standardize_residual_MSN,col="deeppink3")

VIC COMPANY

#Fit the garch model to VIC
garchspec_VIC <- ugarchspec(
  mean.model=list(armaOrder=c(order_VIC[1],order_VIC[3])),
  variance.model=list(model="gjrGARCH",garchOrder = c(1, 1)),
  distribution.model = "sstd")
garchfit_VIC<-ugarchfit(data=VIC_log_diff,spec=garchspec_VIC)
coef(garchfit_VIC)

##            mu           ma1         omega        alpha1         beta1 
##  5.470877e-04 -5.990103e-02  4.178418e-05  3.111150e-01  6.513012e-01 
##        gamma1          skew         shape 
##  7.150895e-02  1.036801e+00  2.881657e+00

#Standardized residuals for VIC
standardize_residual_VIC<-residuals(garchfit_VIC)/sigma(garchfit_VIC)

# Plotting
layout(matrix(c(1,2,3,4),2,2,byrow=TRUE))
acf(standardize_residual_VIC,main="ACF (RESIDUAL-VIC)")
pacf(standardize_residual_VIC,main="PACF (RESIDUAL-VIC)")

plot(standardize_residual_VIC,col="red",main="RESIDUAL-VIC",lwd=1)
qqnorm(standardize_residual_VIC,main="QQ-PLOT VIC")
qqline(standardize_residual_VIC,col="red")

VNM COMPANY

#Fit the garch model to VNM
garchspec_VNM <- ugarchspec(
  mean.model=list(armaOrder=c(order_VNM[1],order_VNM[3])),
  variance.model=list(model="gjrGARCH",garchOrder = c(1, 1)),
  distribution.model = "sstd")
garchfit_VNM<-ugarchfit(data=VNM_log_diff,spec=garchspec_VNM)
coef(garchfit_VNM)

##            mu           ar1           ar2           ar3           ma1 
##  5.276012e-04 -7.635466e-01  8.121380e-01  9.269353e-01  7.622612e-01 
##           ma2           ma3         omega        alpha1         beta1 
## -8.138018e-01 -9.123864e-01  1.537412e-05  1.827079e-01  7.502773e-01 
##        gamma1          skew         shape 
##  3.507683e-02  1.089132e+00  4.686577e+00

#Standardized residuals for MSN
standardize_residual_VNM<-residuals(garchfit_VNM)/sigma(garchfit_VNM)

# Plotting
layout(matrix(c(1,2,3,4),2,2,byrow=TRUE))
acf(standardize_residual_VNM,main="ACF (RESIDUAL-VNM)")
pacf(standardize_residual_VNM,main="PACF (RESIDUAL-VNM)")

plot(standardize_residual_VNM,col="blue",main="RESIDUAL-VNM",lwd=1)
qqnorm(standardize_residual_VNM,main="QQ-PLOT VNM")
qqline(standardize_residual_VNM,col="blue")

# Histogram of each standardized residuals
layout(matrix(c(1,2,3,4),2,2,byrow=TRUE))
hist(standardize_residual_CTG,main='CTG standardized residual',breaks=70,col="darkorange2")
hist(standardize_residual_MSN,main='MSN standardized residual',breaks=70,col="deeppink3")
hist(standardize_residual_VIC,main='VIC standardized residual',breaks=70,col="red")
hist(standardize_residual_VNM,main='VNM standardized residual',breaks=70,col="blue")

JARQUE-BERA TEST

# JARQUE-BERA TEST
## Jarque-Bera test for log-return 
JBtest_CTG_return<-jarque.bera.test(CTG_log_diff)
JBtest_MSN_return<-jarque.bera.test(MSN_log_diff)
JBtest_VIC_return<-jarque.bera.test(VIC_log_diff)
JBtest_VNM_return<-jarque.bera.test(VNM_log_diff)

JBtest_return<-matrix(c(JBtest_CTG_return$statistic,JBtest_CTG_return$p.value,JBtest_MSN_return$statistic,JBtest_MSN_return$p.value,JBtest_VIC_return$statistic,JBtest_VIC_return$p.value,JBtest_VNM_return$statistic,JBtest_VNM_return$p.value),2,4)
dimnames(JBtest_return)<-list(c("Jarque-Bera TS","Jarque-Bera p-value"),c("CTG","MSN","VIC","VNM"))
JBtest_return

##                          CTG      MSN      VIC      VNM
## Jarque-Bera TS      584.6375 390.0399 878.7025 431.3966
## Jarque-Bera p-value   0.0000   0.0000   0.0000   0.0000

LJUNG-BOX TEST

# LJUNG BOX TEST
## LJung- Box test for both log-return and standardized residuals
LJtest_CTG_return<-Box.test(CTG_log_diff)
LJtest_MSN_return<-Box.test(MSN_log_diff)
LJtest_VIC_return<-Box.test(VIC_log_diff)
LJtest_VNM_return<-Box.test(VNM_log_diff)

LJtest_return<-matrix(c(LJtest_CTG_return$statistic,LJtest_CTG_return$p.value,LJtest_MSN_return$statistic,LJtest_MSN_return$p.value,LJtest_VIC_return$statistic,LJtest_VIC_return$p.value,LJtest_VNM_return$statistic,LJtest_VNM_return$p.value),2,4)
dimnames(LJtest_return)<-list(c("LJung-Box TS","LJung-Box p-value"),c("CTG","MSN","VIC","VNM"))
LJtest_return

##                         CTG        MSN       VIC         VNM
## LJung-Box TS      1.4587372 2.86511082 0.1932507 0.008379789
## LJung-Box p-value 0.2271311 0.09051948 0.6602248 0.927062583

LJtest_CTG_residual<-Box.test(standardize_residual_CTG)
LJtest_MSN_residual<-Box.test(standardize_residual_MSN)
LJtest_VIC_residual<-Box.test(standardize_residual_VIC)
LJtest_VNM_residual<-Box.test(standardize_residual_VNM)

LJtest_residual<-matrix(c(LJtest_CTG_residual$statistic,LJtest_CTG_residual$p.value,LJtest_MSN_residual$statistic,LJtest_MSN_residual$p.value,LJtest_VIC_residual$statistic,LJtest_VIC_residual$p.value,LJtest_VNM_residual$statistic,LJtest_VNM_residual$p.value),2,4)
dimnames(LJtest_residual)<-list(c("LJung-Box TS","LJung-Box p-value"),c("CTG","MSN","VIC","VNM"))
LJtest_residual

##                         CTG        MSN        VIC        VNM
## LJung-Box TS      1.3367659 4.01601135 3.00342265 4.05407232
## LJung-Box p-value 0.2476051 0.04507018 0.08308882 0.04406492

ENGLE’S ARCH TEST

# ENGLE'S ARCH TEST
## Engle's test for both log-return and standardized residuals
ARCHtest_CTG_return<-archtest(as.vector(CTG_log_diff))
ARCHtest_MSN_return<-archtest(as.vector(MSN_log_diff))
ARCHtest_VIC_return<-archtest(as.vector(VIC_log_diff))
ARCHtest_VNM_return<-archtest(as.vector(VNM_log_diff))

ARCHtest_return<-matrix(c(ARCHtest_CTG_return$statistic,ARCHtest_CTG_return$p.value,ARCHtest_MSN_return$statistic,ARCHtest_MSN_return$p.value,ARCHtest_VIC_return$statistic,ARCHtest_VIC_return$p.value,ARCHtest_VNM_return$statistic,ARCHtest_VNM_return$p.value),2,4)
dimnames(ARCHtest_return)<-list(c("Engle-ARCH TS","Engle-ARCH p-value"),c("CTG","MSN","VIC","VNM"))
ARCHtest_return

##                             CTG          MSN          VIC          VNM
## Engle-ARCH TS      9.567215e+01 6.471995e+01 4.450997e+01 6.398498e+01
## Engle-ARCH p-value 1.355744e-22 8.633626e-16 2.530656e-11 1.253713e-15

ARCHtest_CTG_residual<-archtest(as.vector(standardize_residual_CTG))
ARCHtest_MSN_residual<-archtest(as.vector(standardize_residual_MSN))
ARCHtest_VIC_residual<-archtest(as.vector(standardize_residual_VIC))
ARCHtest_VNM_residual<-archtest(as.vector(standardize_residual_VNM))

ARCHtest_residual<-matrix(c(ARCHtest_CTG_residual$statistic,ARCHtest_CTG_residual$p.value,ARCHtest_MSN_residual$statistic,ARCHtest_MSN_residual$p.value,ARCHtest_VIC_residual$statistic,ARCHtest_VIC_residual$p.value,ARCHtest_VNM_residual$statistic,ARCHtest_VNM_residual$p.value),2,4)
dimnames(ARCHtest_residual)<-list(c("Engle-ARCH TS","Engle-ARCH p-value"),c("CTG","MSN","VIC","VNM"))

ARCHtest_residual

##                          CTG       MSN       VIC       VNM
## Engle-ARCH TS      0.2397421 0.7934553 0.8945519 0.0115108
## Engle-ARCH p-value 0.6243925 0.3730573 0.3442468 0.9145602

summary_testing<-rbind(rep(1,4),rep(0,4),c(0,1,0,0),rep(1,4),rep(0,4))
dimnames(summary_testing)<-list(c("Jarque-Bera test","LJung-Box test (return)","LJung-Box test (residual)","Engle-ARCH test (return)","Engle-ARCH test (residual)"),c("CTG","MSN","VIC","VNM"))
summary_testing

##                            CTG MSN VIC VNM
## Jarque-Bera test             1   1   1   1
## LJung-Box test (return)      0   0   0   0
## LJung-Box test (residual)    0   1   0   0
## Engle-ARCH test (return)     1   1   1   1
## Engle-ARCH test (residual)   0   0   0   0

PDF FIT OF THE DATA (GENERALIZED PARETO FOR THE TAIL AND NORMAL FOR THE INTERIOR)

x_CTG<-coredata(standardize_residual_CTG)
x_MSN<-coredata(standardize_residual_MSN)
x_VIC<-coredata(standardize_residual_VIC)
x_VNM<-coredata(standardize_residual_VNM)

CTG COMPANY

fit_CTG<-GNG_fit(standardize_residual_CTG, start = c(break1 = -2, break2 =1.5, mean = 0, sd =1,shape1 = 0.1, shape2= 0.1))
fit_CTG

## Fitted composite GPD-Normal-GPD distribution: 
## 
## Breakpoints: -2.453699 1.39936 
## Weights: 0.008083 0.916282 0.075635 
## 
## Parameters: 
##      loc1    scale1    shape1      mean        sd      loc2    scale2    shape2 
##  2.453699  1.429695 -0.508695 -0.055311  0.806488  1.399360  0.817373  0.041094 
## 
## Log-likelihood: -2368.322,  Average log-likelihood: -1.3674

plot(fit_CTG)

MSN COMPANY

fit_MSN<-GNG_fit(standardize_residual_MSN, start = c(break1 = -2, break2 =2, mean = 0, sd = 1,shape1 = 0.2, shape2 = 0.1))
fit_MSN

## Fitted composite GPD-Normal-GPD distribution: 
## 
## Breakpoints: -1.439152 1.250101 
## Weights: 0.046767 0.874711 0.078522 
## 
## Parameters: 
##      loc1    scale1    shape1      mean        sd      loc2    scale2    shape2 
##  1.439152  0.878143 -0.200248 -0.016568  0.657213  1.250101  0.907456 -0.120298 
## 
## Log-likelihood: -2220.541,  Average log-likelihood: -1.2821

plot(fit_MSN)

VIC COMPANY

fit_VIC<-GNG_fit(standardize_residual_VIC, start = c(break1 = -2, break2 =1.5, mean = 0, sd = 1,shape1 = 0.15, shape2 = 0.2))
fit_VIC

## Fitted composite GPD-Normal-GPD distribution: 
## 
## Breakpoints: -1.248442 1.046547 
## Weights: 0.056582 0.854503 0.088915 
## 
## Parameters: 
##      loc1    scale1    shape1      mean        sd      loc2    scale2    shape2 
##  1.248442  0.786275  0.053879 -0.035775  0.585720  1.046547  0.798963 -0.059761 
## 
## Log-likelihood: -2099.176,  Average log-likelihood: -1.212

plot(fit_VIC)

VNM COMPANY

fit_VNM<-GNG_fit(standardize_residual_VNM, start = c(break1 = -2, break2 =2, mean = 0 , sd = 1 ,shape1 = 1, shape2 = 0.5))
fit_VNM

## Fitted composite GPD-Normal-GPD distribution: 
## 
## Breakpoints: -2.486537 1.471217 
## Weights: 0.006928 0.926097 0.066975 
## 
## Parameters: 
##      loc1    scale1    shape1      mean        sd      loc2    scale2    shape2 
##  2.486537  1.155797 -0.299684 -0.066144  0.820428  1.471217  0.833116 -0.106343 
## 
## Log-likelihood: -2364.98,  Average log-likelihood: -1.3655

plot(fit_VNM)

DEFINE THE PDF AND CDF AND INVERSE CDF

# CTG
PDF_CTG<- function(x){d(distribution(fit_CTG),x)}
CDF_CTG<- function(x){p(distribution(fit_CTG),x)}
inverse_CDF_CTG<-function(x){q(distribution(fit_CTG),x)}

# MSN
PDF_MSN<- function(x){d(distribution(fit_MSN),x)}
CDF_MSN<- function(x){p(distribution(fit_MSN),x)}
inverse_CDF_MSN<-function(x){q(distribution(fit_MSN),x)}

# VIC
PDF_VIC<- function(x){d(distribution(fit_VIC),x)}
CDF_VIC<- function(x){p(distribution(fit_VIC),x)}
inverse_CDF_VIC<-function(x){q(distribution(fit_VIC),x)}

# VNM
PDF_VNM<- function(x){d(distribution(fit_VNM),x)}
CDF_VNM<- function(x){p(distribution(fit_VNM),x)}
inverse_CDF_VNM<-function(x){q(distribution(fit_VNM),x)}

# Plot all the graph

# DENSITY DISTRIBUTION FUNCTION
layout(matrix(c(1,2,3,4),2,2,byrow = TRUE))
curve(PDF_CTG,-3,3,col="darkorange2")
curve(PDF_MSN,-3,3,col="deeppink3")
curve(PDF_VIC,-3,3,col="red")
curve(PDF_VNM,-3,3,col="blue")

# CUMMULATIVE DISTRIBUTION FUNCTION
layout(matrix(c(1,2,3,4),2,2,byrow = TRUE))
curve(CDF_CTG,-3,3,col="darkorange2")
curve(CDF_MSN,-3,3,col="deeppink3")
curve(CDF_VIC,-3,3,col="red")
curve(CDF_VNM,-3,3,col="blue")

# INVERSE CDF FUNCTION
layout(matrix(c(1,2,3,4),2,2,byrow = TRUE))
curve(inverse_CDF_CTG,col="darkorange2")
curve(inverse_CDF_MSN,col="deeppink3")
curve(inverse_CDF_VIC,col="red")
curve(inverse_CDF_VNM,col="blue")

TRANSFORM STANDARDIZED RESIDUALS TO THE COPULA SCALE (U DATA)

U_CTG<-CDF_CTG(standardize_residual_CTG)
U_MSN<-CDF_MSN(standardize_residual_MSN)
U_VIC<-CDF_VIC(standardize_residual_VIC)
U_VNM<-CDF_VNM(standardize_residual_VNM)

data<-cbind(U_CTG,U_MSN,U_VIC,U_VNM)

layout(matrix(c(1,2,3,4),2,2,byrow=TRUE))
plot(x_CTG,x_MSN)
plot(x_CTG,x_VIC)
plot(x_CTG,x_VNM)
plot(x_MSN,x_VIC)

plot(x_MSN,x_VNM)
plot(x_VIC,x_VNM)

plot(U_CTG,U_MSN)
plot(U_CTG,U_VIC)

plot(U_CTG,U_VNM)
plot(U_MSN,U_VIC)
plot(U_MSN,U_VNM)
plot(U_VIC,U_VNM)

DEPENDENCE MEASURE BY SPEARMAN RHO AND KENDALL TAU

# Find the correlation between Apple and Facebook by Spearman's rho
rho   <- cor(data, method = "spearman")
rho 

##           U_CTG     U_MSN     U_VIC     U_VNM
## U_CTG 1.0000000 0.2103717 0.2310469 0.2465593
## U_MSN 0.2103717 1.0000000 0.2444499 0.2193768
## U_VIC 0.2310469 0.2444499 1.0000000 0.2054057
## U_VNM 0.2465593 0.2193768 0.2054057 1.0000000

# Find the correlation between Apple and Facebook by Kendall tau
kendall<-cor(data, method = "kendall")
kendall

##           U_CTG     U_MSN     U_VIC     U_VNM
## U_CTG 1.0000000 0.1414446 0.1566663 0.1658348
## U_MSN 0.1414446 1.0000000 0.1657134 0.1484544
## U_VIC 0.1566663 0.1657134 1.0000000 0.1383667
## U_VNM 0.1658348 0.1484544 0.1383667 1.0000000

FIT THE COPULA TO THE DATA AND FIND COPULA PARAMETERS

Gaussian Copula

# we fit all the Copula by Maximum likelihood method (ml)
X<-as.matrix(data)
fit_Gaussian <- fitCopula(normalCopula(dim=4), X, method = 'ml')
summary(fit_Gaussian)

## Call: fitCopula(copula, data = data, method = "ml")
## Fit based on "maximum likelihood" and 1732 4-dimensional observations.
## Normal copula, dim. d = 4 
##       Estimate Std. Error
## rho.1   0.2461      0.012
## The maximized loglikelihood is 242 
## Optimization converged
## Number of loglikelihood evaluations:
## function gradient 
##        6        6

coef(fit_Gaussian)

##     rho.1 
## 0.2461032

t-Student Copula

fit_tStudent <- fitCopula(tCopula(dim=4),X, method = 'ml')
summary(fit_tStudent)

## Call: fitCopula(copula, data = data, method = "ml")
## Fit based on "maximum likelihood" and 1732 4-dimensional observations.
## t-copula, dim. d = 4 
##       Estimate Std. Error
## rho.1   0.2451      0.013
## df     28.2030      8.770
## The maximized loglikelihood is 248 
## Optimization converged
## Number of loglikelihood evaluations:
## function gradient 
##       25       25

coef(fit_tStudent)

##      rho.1         df 
##  0.2451165 28.2030387

Clayton Copula

fit_Clayton <- fitCopula(claytonCopula(dim=4), X, method = 'itau')
summary(fit_Clayton)

## Call: fitCopula(copula, data = data, method = "itau")
## Fit based on "inversion of Kendall's tau" and 1732 4-dimensional observations.
## Clayton copula, dim. d = 4 
##       Estimate Std. Error
## alpha    0.361      0.026

coef(fit_Clayton)

##     alpha 
## 0.3609571

Gumbel Copula

fit_Gumbel <- fitCopula(gumbelCopula(dim=4), X, method = 'itau')
summary(fit_Gumbel)

## Call: fitCopula(copula, data = data, method = "itau")
## Fit based on "inversion of Kendall's tau" and 1732 4-dimensional observations.
## Gumbel copula, dim. d = 4 
##       Estimate Std. Error
## alpha     1.18      0.013

coef(fit_Gumbel)

##    alpha 
## 1.180479

Frank Copula

fit_Frank <- fitCopula(frankCopula(dim=4), X, method = 'itau')
summary(fit_Frank)

## Call: fitCopula(copula, data = data, method = "itau")
## Fit based on "inversion of Kendall's tau" and 1732 4-dimensional observations.
## Frank copula, dim. d = 4 
##       Estimate Std. Error
## alpha    1.402      0.005

coef(fit_Frank)

##    alpha 
## 1.401783

CONSTRUCT THE COPULA

## Note: 
## Sigma of Gaussian Copula and T-Student Copula should be around 0.205-0.2465 
## Tau value of Clayton, Gumbel and Frank Copula should be around 0.138-0.165

# Gaussian Copula
Gaussian_model<-normalCopula(coef(fit_Gaussian) ,dim=4)
getSigma(Gaussian_model)

##           [,1]      [,2]      [,3]      [,4]
## [1,] 1.0000000 0.2461032 0.2461032 0.2461032
## [2,] 0.2461032 1.0000000 0.2461032 0.2461032
## [3,] 0.2461032 0.2461032 1.0000000 0.2461032
## [4,] 0.2461032 0.2461032 0.2461032 1.0000000

# Student T Copula
T_model<-tCopula(coef(fit_tStudent)[1] ,dim=4,df=coef(fit_tStudent)[2])
getSigma(T_model)

##           [,1]      [,2]      [,3]      [,4]
## [1,] 1.0000000 0.2451165 0.2451165 0.2451165
## [2,] 0.2451165 1.0000000 0.2451165 0.2451165
## [3,] 0.2451165 0.2451165 1.0000000 0.2451165
## [4,] 0.2451165 0.2451165 0.2451165 1.0000000

# Clayton Copula
Clayton_model<-claytonCopula(coef(fit_Clayton),dim=4)
tau(Clayton_model)

##     alpha 
## 0.1528859

# Gumbel Copula
Gumbel_model<-gumbelCopula( coef(fit_Gumbel),dim=4)
tau(Gumbel_model)

##     alpha 
## 0.1528859

# Frank Copula
Frank_model<-frankCopula(coef(fit_Frank),dim=4)
tau(Frank_model)

##     alpha 
## 0.1527916

SIMULATION OF EACH TYPE OF COPULA

n=1000
m=nrow(data)

# Gaussian Copula
CTG_matrix_G<-matrix(rep(0,n*m),nrow=m)
MSN_matrix_G<-matrix(rep(0,n*m),nrow=m)
VIC_matrix_G<-matrix(rep(0,n*m),nrow=m)
VNM_matrix_G<-matrix(rep(0,n*m),nrow=m)

# T-Student Copula
CTG_matrix_T<-matrix(rep(0,n*m),nrow=m)
MSN_matrix_T<-matrix(rep(0,n*m),nrow=m)
VIC_matrix_T<-matrix(rep(0,n*m),nrow=m)
VNM_matrix_T<-matrix(rep(0,n*m),nrow=m)

# Clayton Copula
CTG_matrix_C<-matrix(rep(0,n*m),nrow=m)
MSN_matrix_C<-matrix(rep(0,n*m),nrow=m)
VIC_matrix_C<-matrix(rep(0,n*m),nrow=m)
VNM_matrix_C<-matrix(rep(0,n*m),nrow=m)

# Gumbel Copula
CTG_matrix_Gum<-matrix(rep(0,n*m),nrow=m)
MSN_matrix_Gum<-matrix(rep(0,n*m),nrow=m)
VIC_matrix_Gum<-matrix(rep(0,n*m),nrow=m)
VNM_matrix_Gum<-matrix(rep(0,n*m),nrow=m)

# Frank Copula
CTG_matrix_F<-matrix(rep(0,n*m),nrow=m)
MSN_matrix_F<-matrix(rep(0,n*m),nrow=m)
VIC_matrix_F<-matrix(rep(0,n*m),nrow=m)
VNM_matrix_F<-matrix(rep(0,n*m),nrow=m)

GAUSSIAN COPULA

for(i in 1:n){
    set.seed(i)
    U<- rCopula(m, copula = Gaussian_model)
    CTG_matrix_G[,i]<-inverse_CDF_CTG(U[,1])
    MSN_matrix_G[,i]<-inverse_CDF_MSN(U[,2])
    VIC_matrix_G[,i]<-inverse_CDF_VIC(U[,3])
    VNM_matrix_G[,i]<-inverse_CDF_VNM(U[,4])
}

# Turn to xts objects
SR_CTG_G<-xts(x =CTG_matrix_G,order.by=index(standardize_residual_CTG))
SR_MSN_G<-xts(x =MSN_matrix_G,order.by=index(standardize_residual_MSN))
SR_VIC_G<-xts(x =VIC_matrix_G,order.by=index(standardize_residual_VIC))
SR_VNM_G<-xts(x =VNM_matrix_G,order.by=index(standardize_residual_VNM))

# Plotting
layout(matrix(c(1,2,1,2),2,2))
plot(SR_CTG_G,ylim=c(-7,7),main="Simulated S.RESIDUAL-CTG (Gaussian)",lwd=0.5)
plot(SR_MSN_G,ylim=c(-7,7),main="Simulated S.RESIDUAL-MSN (Gaussian)",lwd=0.5)

plot(SR_VIC_G,ylim=c(-7,7),main="Simulated S.RESIDUAL-VIC (Gaussian)",lwd=0.5)
plot(SR_VNM_G,ylim=c(-7,7),main="Simulated S.RESIDUAL-VNM (Gaussian)",lwd=0.5)

T-STUDENT COPULA

for(i in 1:n){
    set.seed(i)
    U<- rCopula(m, copula = T_model)
    CTG_matrix_T[,i]<-inverse_CDF_CTG(U[,1])
    MSN_matrix_T[,i]<-inverse_CDF_MSN(U[,2])
    VIC_matrix_T[,i]<-inverse_CDF_VIC(U[,3])
    VNM_matrix_T[,i]<-inverse_CDF_VNM(U[,4])
}

# Turn to xts objects
SR_CTG_T<-xts(x =CTG_matrix_T,order.by=index(standardize_residual_CTG))
SR_MSN_T<-xts(x =MSN_matrix_T,order.by=index(standardize_residual_MSN))
SR_VIC_T<-xts(x =VIC_matrix_T,order.by=index(standardize_residual_VIC))
SR_VNM_T<-xts(x =VNM_matrix_T,order.by=index(standardize_residual_VNM))

# Plotting
layout(matrix(c(1,2,1,2),2,2))
plot(SR_CTG_T,ylim=c(-7,7),main="Simulated S.RESIDUAL-CTG (T-Student)",lwd=0.5)
plot(SR_MSN_T,ylim=c(-7,7),main="Simulated S.RESIDUAL-MSN (T-Student)",lwd=0.5)

plot(SR_VIC_T,ylim=c(-7,7),main="Simulated S.RESIDUAL-VIC (T-Student)",lwd=0.5)
plot(SR_VNM_T,ylim=c(-7,7),main="Simulated S.RESIDUAL-VNM (T-Student)",lwd=0.5)

CLAYTON COPULA

for(i in 1:n){
  set.seed(i)
    U<- rCopula(m, copula = Clayton_model)
    CTG_matrix_C[,i]<-inverse_CDF_CTG(U[,1])
    MSN_matrix_C[,i]<-inverse_CDF_MSN(U[,2])
    VIC_matrix_C[,i]<-inverse_CDF_VIC(U[,3])
    VNM_matrix_C[,i]<-inverse_CDF_VNM(U[,4])
}

# Turn to xts objects
SR_CTG_C<-xts(x =CTG_matrix_C,order.by=index(standardize_residual_CTG))
SR_MSN_C<-xts(x =MSN_matrix_C,order.by=index(standardize_residual_MSN))
SR_VIC_C<-xts(x =VIC_matrix_C,order.by=index(standardize_residual_VIC))
SR_VNM_C<-xts(x =VNM_matrix_C,order.by=index(standardize_residual_VNM))

# Plotting
layout(matrix(c(1,2,1,2),2,2))
plot(SR_CTG_C,ylim=c(-7,7),main="Simulated S.RESIDUAL-CTG (Clayton)",lwd=0.5)
plot(SR_MSN_C,ylim=c(-7,7),main="Simulated S.RESIDUAL-MSN (Clayton)",lwd=0.5)

plot(SR_VIC_C,ylim=c(-7,7),main="Simulated S.RESIDUAL-VIC (Clayton)",lwd=0.5)
plot(SR_VNM_C,ylim=c(-7,7),main="Simulated S.RESIDUAL-VNM (Clayton)",lwd=0.5)

GUMBEL COPULA

for(i in 1:n){
    set.seed(i)
    U<- rCopula(m, copula = Gumbel_model)
    CTG_matrix_Gum[,i]<-inverse_CDF_CTG(U[,1])
    MSN_matrix_Gum[,i]<-inverse_CDF_MSN(U[,2])
    VIC_matrix_Gum[,i]<-inverse_CDF_VIC(U[,3])
    VNM_matrix_Gum[,i]<-inverse_CDF_VNM(U[,4])
}

# Turn to xts objects
SR_CTG_Gum<-xts(x =CTG_matrix_Gum,order.by=index(standardize_residual_CTG))
SR_MSN_Gum<-xts(x =MSN_matrix_Gum,order.by=index(standardize_residual_MSN))
SR_VIC_Gum<-xts(x =VIC_matrix_Gum,order.by=index(standardize_residual_VIC))
SR_VNM_Gum<-xts(x =VNM_matrix_Gum,order.by=index(standardize_residual_VNM))

# Plotting
layout(matrix(c(1,2,1,2),2,2))
plot(SR_CTG_Gum,ylim=c(-7,7),main="Simulated S.RESIDUAL-CTG (Gumbel)",lwd=0.5)
plot(SR_MSN_Gum,ylim=c(-7,7),main="Simulated S.RESIDUAL-MSN (Gumbel)",lwd=0.5)

plot(SR_VIC_Gum,ylim=c(-7,7),main="Simulated S.RESIDUAL-VIC (Gumbel)",lwd=0.5)
plot(SR_VNM_Gum,ylim=c(-7,7),main="Simulated S.RESIDUAL-VNM (Gumbel)",lwd=0.5)

FRANK COPULA

for(i in 1:n){
    set.seed(i)
    U<- rCopula(m, copula = Frank_model)
    CTG_matrix_F[,i]<-inverse_CDF_CTG(U[,1])
    MSN_matrix_F[,i]<-inverse_CDF_MSN(U[,2])
    VIC_matrix_F[,i]<-inverse_CDF_VIC(U[,3])
    VNM_matrix_F[,i]<-inverse_CDF_VNM(U[,4])
}

# Turn to xts objects
SR_CTG_F<-xts(x =CTG_matrix_F,order.by=index(standardize_residual_CTG))
SR_MSN_F<-xts(x =MSN_matrix_F,order.by=index(standardize_residual_MSN))
SR_VIC_F<-xts(x =VIC_matrix_F,order.by=index(standardize_residual_VIC))
SR_VNM_F<-xts(x =VNM_matrix_F,order.by=index(standardize_residual_VNM))

# Plotting
layout(matrix(c(1,2,1,2),2,2))
plot(SR_CTG_F,ylim=c(-7,7),main="Simulated S.RESIDUAL-CTG (Frank)",lwd=0.5)
plot(SR_MSN_F,ylim=c(-7,7),main="Simulated S.RESIDUAL-MSN (Frank)",lwd=0.5)

plot(SR_VIC_F,ylim=c(-7,7),main="Simulated S.RESIDUAL-VIC (Frank)",lwd=0.5)
plot(SR_VNM_F,ylim=c(-7,7),main="Simulated S.RESIDUAL-VNM (Frank)",lwd=0.5)

TURN STANDARDIZED RESIDUAL TO LOG-RETURN (REINTRODUCE ARMA-GARCH MODEL)

# Create a function to turn standardized residual to log-return
return <- function(model,SR ,num_col)
  {
sigma_matrix<-coredata(sigma(model))
diagonal_sigma<-diag(as.numeric(sigma_matrix),nrow=length(sigma_matrix))
simulated_residuals<-diagonal_sigma%*%SR
simulated_log_returns<-simulated_residuals+matrix(rep(as.numeric(coredata(fitted(model))),num_col),ncol=num_col)
simulated_log_returns<-xts(simulated_log_returns,order.by=index(sigma(model)))
simulated_log_returns
}

Gaussian Copula

return_CTG_G<- return(garchfit_CTG,SR_CTG_G,num_col=n)
return_MSN_G<- return(garchfit_MSN,SR_MSN_G,num_col=n)
return_VIC_G<- return(garchfit_VIC,SR_VIC_G,num_col=n)
return_VNM_G<- return(garchfit_VNM,SR_VNM_G,num_col=n)

# Plotting
layout(matrix(c(1,2,1,2),2,2))
plot(return_CTG_G,main="Simulated RETURN-CTG (Gaussian)",lwd=0.5)
plot(return_MSN_G,main="Simulated RETURN-MSN (Gaussian)",lwd=0.5)

plot(return_VIC_G,main="Simulated RETURN-VIC (Gaussian)",lwd=0.5)
plot(return_VNM_G,main="Simulated RETURN-VNM (Gaussian)",lwd=0.5)

T-student Copula

return_CTG_T<- return(garchfit_CTG,SR_CTG_T,num_col=n)
return_MSN_T<- return(garchfit_MSN,SR_MSN_T,num_col=n)
return_VIC_T<- return(garchfit_VIC,SR_VIC_T,num_col=n)
return_VNM_T<- return(garchfit_VNM,SR_VNM_T,num_col=n)

# Plotting
layout(matrix(c(1,2,1,2),2,2))
plot(return_CTG_T,main="Simulated RETURN-CTG (T-Student)",lwd=0.5)
plot(return_MSN_T,main="Simulated RETURN-MSN (T-Student)",lwd=0.5)

plot(return_VIC_T,main="Simulated RETURN-VIC (T-Student)",lwd=0.5)
plot(return_VNM_T,main="Simulated RETURN-VNM (T-Student)",lwd=0.5)

Clayton Copula

return_CTG_C<- return(garchfit_CTG,SR_CTG_C,num_col=n)
return_MSN_C<- return(garchfit_MSN,SR_MSN_C,num_col=n)
return_VIC_C<- return(garchfit_VIC,SR_VIC_C,num_col=n)
return_VNM_C<- return(garchfit_VNM,SR_VNM_C,num_col=n)

# Plotting
layout(matrix(c(1,2,1,2),2,2))
plot(return_CTG_C,main="Simulated RETURN-CTG (Clayton)",lwd=0.5)
plot(return_MSN_C,main="Simulated RETURN-MSN (Clayton)",lwd=0.5)

plot(return_VIC_C,main="Simulated RETURN-VIC (Clayton)",lwd=0.5)
plot(return_VNM_C,main="Simulated RETURN-VNM (Clayton)",lwd=0.5)

Gumbel Copula

return_CTG_Gum<- return(garchfit_CTG,SR_CTG_Gum,num_col=n)
return_MSN_Gum<- return(garchfit_MSN,SR_MSN_Gum,num_col=n)
return_VIC_Gum<- return(garchfit_VIC,SR_VIC_Gum,num_col=n)
return_VNM_Gum<- return(garchfit_VNM,SR_VNM_Gum,num_col=n)

# Plotting
layout(matrix(c(1,2,1,2),2,2))
plot(return_CTG_Gum,main="Simulated RETURN-CTG (Gumbel)",lwd=0.5)
plot(return_MSN_Gum,main="Simulated RETURN-MSN (Gumbel)",lwd=0.5)

plot(return_VIC_Gum,main="Simulated RETURN-VIC (Gumbel)",lwd=0.5)
plot(return_VNM_Gum,main="Simulated RETURN-VNM (Gumbel)",lwd=0.5)

Frank Copula

return_CTG_F<- return(garchfit_CTG,SR_CTG_F,num_col=n)
return_MSN_F<- return(garchfit_MSN,SR_MSN_F,num_col=n)
return_VIC_F<- return(garchfit_VIC,SR_VIC_F,num_col=n)
return_VNM_F<- return(garchfit_VNM,SR_VNM_F,num_col=n)

# Plotting
layout(matrix(c(1,2,1,2),2,2))
plot(return_CTG_F,main="Simulated RETURN-CTG (Frank)",lwd=0.5)
plot(return_MSN_F,main="Simulated RETURN-MSN (Frank)",lwd=0.5)

plot(return_VIC_F,main="Simulated RETURN-VIC (Frank)",lwd=0.5)
plot(return_VNM_F,main="Simulated RETURN-VNM (Frank)",lwd=0.5)

VAR CALCULAION

weight<-c(0.25,0.25,0.25,0.25)

real_port_return<-(weight[1]*CTG_log_diff)+(weight[2]*MSN_log_diff)+(weight[3]*VIC_log_diff)+(weight[4]*VNM_log_diff)

# Plotting
plot(real_port_return,ylim=c(-0.04,0.04),main="real portfolio return",lwd=1,col="darkblue")

Gaussian Copula

# Calculate simulated portfolio return
port_return_G<-weight[1]*return_CTG_G +weight[2]*return_MSN_G +weight[3]*return_VIC_G +weight[4]*return_VNM_G

# VaR 99
VaR99_return_G<-xts(matrix(rep(1,m),ncol=1),order.by=index((port_return_G)))
for(i in 1:m){
VaR99_return_G[i,]<-quantile(port_return_G[i,],p=0.01)
}

# VaR 95
VaR95_return_G<-xts(matrix(rep(1,m),ncol=1),order.by=index((port_return_G)))
for(i in 1:m){
VaR95_return_G[i,]<-quantile(port_return_G[i,],p=0.05)
}

# VaR 90
VaR90_return_G<-xts(matrix(rep(1,m),ncol=1),order.by=index((port_return_G)))
for(i in 1:m){
VaR90_return_G[i,]<-quantile(port_return_G[i,],p=0.1)
}

# Plotting
layout(matrix(c(1,2,1,2),2,2))
plot(port_return_G,ylim=c(-0.1,0.15),main="Simulated port_return (Gaussian)",lwd=0.5)
VaRplot(0.01,real_port_return,VaR99_return_G)

VaRplot(0.05,real_port_return,VaR95_return_G)
VaRplot(0.1,real_port_return,VaR90_return_G)

layout(matrix(c(1),1,1))
plot(cbind(real_port_return,VaR99_return_G,VaR95_return_G,VaR90_return_G),col=c("black","deeppink3","darkorange2","blue"),main="GAUSSIAN COPULA",lwd=0.7)

T_Student Copula

# Calculate simulated portfolio return
port_return_T<-weight[1]*return_CTG_T +weight[2]*return_MSN_T +weight[3]*return_VIC_T +weight[4]*return_VNM_T

# VaR 99
VaR99_return_T<-xts(matrix(rep(1,m),ncol=1),order.by=index((port_return_T)))
for(i in 1:m){
VaR99_return_T[i,]<-quantile(port_return_T[i,],p=0.01)
}

# VaR 95
VaR95_return_T<-xts(matrix(rep(1,m),ncol=1),order.by=index((port_return_T)))
for(i in 1:m){
VaR95_return_T[i,]<-quantile(port_return_T[i,],p=0.05)
}

# VaR 90
VaR90_return_T<-xts(matrix(rep(1,m),ncol=1),order.by=index((port_return_T)))
for(i in 1:m){
VaR90_return_T[i,]<-quantile(port_return_T[i,],p=0.1)
}

# Plotting
layout(matrix(c(1,2,1,2),2,2))
plot(port_return_T,ylim=c(-0.1,0.1),main="Simulated port_return (T-Student)",lwd=0.5)
VaRplot(0.01,real_port_return,VaR99_return_T)

VaRplot(0.05,real_port_return,VaR95_return_T)
VaRplot(0.1,real_port_return,VaR90_return_T)

layout(matrix(c(1),1,1))
plot(cbind(real_port_return,VaR99_return_T,VaR95_return_T,VaR90_return_T),col=c("black","deeppink3","darkorange2","blue"),main="T-STUDENT COPULA",lwd=0.7)

Clayton Copula

# Calculate simulated portfolio return
port_return_C<-weight[1]*return_CTG_C +weight[2]*return_MSN_C +weight[3]*return_VIC_C +weight[4]*return_VNM_C

# VaR 99
VaR99_return_C<-xts(matrix(rep(1,m),ncol=1),order.by=index((port_return_C)))
for(i in 1:m){
VaR99_return_C[i,]<-quantile(port_return_C[i,],p=0.01)
}

# VaR 95
VaR95_return_C<-xts(matrix(rep(1,m),ncol=1),order.by=index((port_return_C)))
for(i in 1:m){
VaR95_return_C[i,]<-quantile(port_return_C[i,],p=0.05)
}

# VaR 90
VaR90_return_C<-xts(matrix(rep(1,m),ncol=1),order.by=index((port_return_C)))
for(i in 1:m){
VaR90_return_C[i,]<-quantile(port_return_C[i,],p=0.1)
}

# Plotting
layout(matrix(c(1,2,1,2),2,2))
plot(port_return_C,ylim=c(-0.13,0.1),main="Simulated port_return (Clayton)",lwd=0.5)
VaRplot(0.01,real_port_return,VaR99_return_C)

VaRplot(0.05,real_port_return,VaR95_return_C)
VaRplot(0.1,real_port_return,VaR90_return_C)

layout(matrix(c(1),1,1))
plot(cbind(real_port_return,VaR99_return_C,VaR95_return_C,VaR90_return_C),col=c("black","deeppink3","darkorange2","blue"),main="CLAYTON COPULA",lwd=0.7)

Gumbel Copula

# Calculate simulated portfolio return
port_return_Gum<-weight[1]*return_CTG_Gum +weight[2]*return_MSN_Gum +weight[3]*return_VIC_Gum +weight[4]*return_VNM_Gum

# VaR 99
VaR99_return_Gum<-xts(matrix(rep(1,m),ncol=1),order.by=index((port_return_Gum)))
for(i in 1:m){
VaR99_return_Gum[i,]<-quantile(port_return_Gum[i,],p=0.01)
}

# VaR 95
VaR95_return_Gum<-xts(matrix(rep(1,m),ncol=1),order.by=index((port_return_Gum)))
for(i in 1:m){
VaR95_return_Gum[i,]<-quantile(port_return_Gum[i,],p=0.05)
}

# VaR 90
VaR90_return_Gum<-xts(matrix(rep(1,m),ncol=1),order.by=index((port_return_Gum)))
for(i in 1:m){
VaR90_return_Gum[i,]<-quantile(port_return_Gum[i,],p=0.1)
}

# Plotting
layout(matrix(c(1,2,1,2),2,2))
plot(port_return_Gum,ylim=c(-0.1,0.17),main="Simulated port_return (Gumbel)",lwd=0.5)
VaRplot(0.01,real_port_return,VaR99_return_Gum)

VaRplot(0.05,real_port_return,VaR95_return_Gum)
VaRplot(0.1,real_port_return,VaR90_return_Gum)

layout(matrix(c(1),1,1))
plot(cbind(real_port_return,VaR99_return_Gum,VaR95_return_Gum,VaR90_return_Gum),col=c("black","deeppink3","darkorange2","blue"),main="GUMBEL COPULA",lwd=0.7)

Frank Copula

# Calculate simulated portfolio return
port_return_F<-weight[1]*return_CTG_F +weight[2]*return_MSN_F +weight[3]*return_VIC_F +weight[4]*return_VNM_F

# VaR 99
VaR99_return_F<-xts(matrix(rep(1,m),ncol=1),order.by=index((port_return_F)))
for(i in 1:m){
VaR99_return_F[i,]<-quantile(port_return_F[i,],p=0.01)
}

# VaR 95
VaR95_return_F<-xts(matrix(rep(1,m),ncol=1),order.by=index((port_return_F)))
for(i in 1:m){
VaR95_return_F[i,]<-quantile(port_return_F[i,],p=0.05)
}

# VaR 90
VaR90_return_F<-xts(matrix(rep(1,m),ncol=1),order.by=index((port_return_F)))
for(i in 1:m){
VaR90_return_F[i,]<-quantile(port_return_F[i,],p=0.1)
}

# Plotting
layout(matrix(c(1,2,1,2),2,2))
plot(port_return_F,ylim=c(-0.08,0.1),main="Simulated port_return (Frank)",lwd=0.5)
VaRplot(0.01,real_port_return,VaR99_return_F)

VaRplot(0.05,real_port_return,VaR95_return_F)
VaRplot(0.1,real_port_return,VaR90_return_F)

layout(matrix(c(1),1,1))
plot(cbind(real_port_return,VaR99_return_F,VaR95_return_F,VaR90_return_F),col=c("black","deeppink3","darkorange2","blue"),main="FRANK COPULA",lwd=0.7)

PLOTTING OF VAR 99%, 95% AND 90% BY EACH TYPE OF COPULA

# Plotting for each vaR 99%, 95% and 90%

## black: real portfolio return 
## blue:  VaR (Gaussuan Copula)
## pink:  VaR (T-Student Copula)
## orange:VaR (Clayton Copula)
## red:   VaR (Gumbel Copula)
## purple:VaR (Frank Copula)

layout(matrix(c(1),1,1))

# VaR 99%
plot(plot(cbind(real_port_return,VaR99_return_G,VaR99_return_T,VaR99_return_C,VaR99_return_Gum,VaR99_return_F),col=c("black","blue","pink","darkorange2","red","darkmagenta"),main="VAR 99%",lwd=0.7))

# VaR 95%
plot(plot(cbind(real_port_return,VaR95_return_G,VaR95_return_T,VaR95_return_C,VaR95_return_Gum,VaR95_return_F),col=c("black","blue","pink","darkorange2","red","darkmagenta"),main="VAR 95%",lwd=0.7))

# VaR 90%
plot(plot(cbind(real_port_return,VaR90_return_G,VaR90_return_T,VaR90_return_C,VaR90_return_Gum,VaR90_return_F),col=c("black","blue","pink","darkorange2","red","darkmagenta"),main="VAR 90%",lwd=0.7))

KUPIEC’S PROPORTION OF FAILURE TEST (UNCONDITIONAL) AND CHRISTOFFERSEN’S INDEPENDENCE TEST (CONDITIONAL)

## we test the hypothesis at anpha= 5% significance level 

# Gaussian Copula
test1_99_G<-VaRTest(0.01,real_port_return,VaR99_return_G)
test1_95_G<-VaRTest(0.05,real_port_return,VaR95_return_G)
test1_90_G<-VaRTest(0.1,real_port_return,VaR90_return_G)

# T Copula
test1_99_T<-VaRTest(0.01,real_port_return,VaR99_return_T)
test1_95_T<-VaRTest(0.05,real_port_return,VaR95_return_T)
test1_90_T<-VaRTest(0.1,real_port_return,VaR90_return_T)

# Clayton Copula
test1_99_C<-VaRTest(0.01,real_port_return,VaR99_return_C)
test1_95_C<-VaRTest(0.05,real_port_return,VaR95_return_C)
test1_90_C<-VaRTest(0.1,real_port_return,VaR90_return_C)

# Gumbell Copula
test1_99_Gum<-VaRTest(0.01,real_port_return,VaR99_return_Gum)
test1_95_Gum<-VaRTest(0.05,real_port_return,VaR95_return_Gum)
test1_90_Gum<-VaRTest(0.1,real_port_return,VaR90_return_Gum)

# Frank Copula
test1_99_F<-VaRTest(0.01,real_port_return,VaR99_return_F)
test1_95_F<-VaRTest(0.05,real_port_return,VaR95_return_F)
test1_90_F<-VaRTest(0.1,real_port_return,VaR90_return_F)

VAR DURATION TEST OF CHRISTOFFERSEN

## we test the hypothesis at anpha= 5% significance level 

# Gaussian Copula
test2_99_G<-VaRDurTest(0.01,real_port_return,VaR99_return_G)
test2_95_G<-VaRDurTest(0.05,real_port_return,VaR95_return_G)
test2_90_G<-VaRDurTest(0.1,real_port_return,VaR90_return_G)

# T Copula
test2_99_T<-VaRDurTest(0.01,real_port_return,VaR99_return_T)
test2_95_T<-VaRDurTest(0.05,real_port_return,VaR95_return_T)
test2_90_T<-VaRDurTest(0.1,real_port_return,VaR90_return_T)

# Clayton Copula
test2_99_C<-VaRDurTest(0.01,real_port_return,VaR99_return_C)
test2_95_C<-VaRDurTest(0.05,real_port_return,VaR95_return_C)
test2_90_C<-VaRDurTest(0.1,real_port_return,VaR90_return_C)

# Gumbell Copula
test2_99_Gum<-VaRDurTest(0.01,real_port_return,VaR99_return_Gum)
test2_95_Gum<-VaRDurTest(0.05,real_port_return,VaR95_return_Gum)
test2_90_Gum<-VaRDurTest(0.1,real_port_return,VaR90_return_Gum)

# Frank Copula
test2_99_F<-VaRDurTest(0.01,real_port_return,VaR99_return_F)
test2_95_F<-VaRDurTest(0.05,real_port_return,VaR95_return_F)
test2_90_F<-VaRDurTest(0.1,real_port_return,VaR90_return_F)

# Gaussian Copula
Kupiec_G<-c(test1_99_G$uc.LRp,test1_99_G$uc.Decision,test1_95_G$uc.LRp,test1_95_G$uc.Decision,test1_90_G$uc.LRp,test1_90_G$uc.Decision)
CI_G<-c(test1_99_G$cc.LRp,test1_99_G$cc.Decision,test1_95_G$cc.LRp,test1_95_G$cc.Decision,test1_90_G$cc.LRp,test1_90_G$cc.Decision)
CD_G<-c(test2_99_G$LRp,test2_99_G$Decision,test2_95_G$LRp,test2_95_G$Decision,test2_90_G$LRp,test2_90_G$Decision)

rbind(Kupiec_G,CI_G,CD_G)

##          [,1]                [,2]                [,3]               
## Kupiec_G "0.127451838037815" "Fail to Reject H0" "0.709472556532043"
## CI_G     "0.223335474037035" "Fail to Reject H0" "0.30928864464668" 
## CD_G     "0.775083080847009" "Fail to Reject H0" "0.577519824048201"
##          [,4]                [,5]                 [,6]               
## Kupiec_G "Fail to Reject H0" "0.161866230228559"  "Fail to Reject H0"
## CI_G     "Fail to Reject H0" "0.0592657871482771" "Fail to Reject H0"
## CD_G     "Fail to Reject H0" "0.834491210799811"  "Fail to Reject H0"

# T-Student Copula
Kupiec_T<-c(test1_99_T$uc.LRp,test1_99_T$uc.Decision,test1_95_T$uc.LRp,test1_95_T$uc.Decision,test1_90_T$uc.LRp,test1_90_T$uc.Decision)
CI_T<-c(test1_99_T$cc.LRp,test1_99_T$cc.Decision,test1_95_T$cc.LRp,test1_95_T$cc.Decision,test1_90_T$cc.LRp,test1_90_T$cc.Decision)
CD_T<-c(test2_99_T$LRp,test2_99_T$Decision,test2_95_T$LRp,test2_95_T$Decision,test2_90_T$LRp,test2_90_T$Decision)

rbind(Kupiec_T,CI_T,CD_T)

##          [,1]                [,2]                [,3]               
## Kupiec_T "0.127451838037815" "Fail to Reject H0" "0.877642134656625"
## CI_T     "0.223335474037035" "Fail to Reject H0" "0.275176713173505"
## CD_T     "0.813686514479523" "Fail to Reject H0" "0.671021081610266"
##          [,4]                [,5]                 [,6]               
## Kupiec_T "Fail to Reject H0" "0.217239021185684"  "Fail to Reject H0"
## CI_T     "Fail to Reject H0" "0.0549677764086901" "Fail to Reject H0"
## CD_T     "Fail to Reject H0" "0.807519557605753"  "Fail to Reject H0"

# Clayton Copula
Kupiec_C<-c(test1_99_C$uc.LRp,test1_99_C$uc.Decision,test1_95_C$uc.LRp,test1_95_C$uc.Decision,test1_90_C$uc.LRp,test1_90_C$uc.Decision)
CI_C<-c(test1_99_C$cc.LRp,test1_99_C$cc.Decision,test1_95_C$cc.LRp,test1_95_C$cc.Decision,test1_90_C$cc.LRp,test1_90_C$cc.Decision)
CD_C<-c(test2_99_C$LRp,test2_99_C$Decision,test2_95_C$LRp,test2_95_C$Decision,test2_90_C$LRp,test2_90_C$Decision)

rbind(Kupiec_C,CI_C,CD_C)

##          [,1]                [,2]                [,3]               
## Kupiec_C "0.275097996415791" "Fail to Reject H0" "0.233068424544467"
## CI_C     "0.49960103611137"  "Fail to Reject H0" "0.458792393527457"
## CD_C     "0.716424388560727" "Fail to Reject H0" "0.776832763324041"
##          [,4]                [,5]                 [,6]               
## Kupiec_C "Fail to Reject H0" "0.0994796964751316" "Fail to Reject H0"
## CI_C     "Fail to Reject H0" "0.113082421436528"  "Fail to Reject H0"
## CD_C     "Fail to Reject H0" "0.990386754241045"  "Fail to Reject H0"

# Gumbel Copula
Kupiec_Gum<-c(test1_99_Gum$uc.LRp,test1_99_Gum$uc.Decision,test1_95_Gum$uc.LRp,test1_95_Gum$uc.Decision,test1_90_Gum$uc.LRp,test1_90_Gum$uc.Decision)
CI_Gum<-c(test1_99_Gum$cc.LRp,test1_99_Gum$cc.Decision,test1_95_Gum$cc.LRp,test1_95_Gum$cc.Decision,test1_90_Gum$cc.LRp,test1_90_Gum$cc.Decision)
CD_Gum<-c(test2_99_Gum$LRp,test2_99_Gum$Decision,test2_95_Gum$LRp,test2_95_Gum$Decision,test2_90_Gum$LRp,test2_90_Gum$Decision)

rbind(Kupiec_Gum,CI_Gum,CD_Gum)

##            [,1]                   [,2]                [,3]                
## Kupiec_Gum "8.1483655644754e-05"  "Reject H0"         "0.0384915877825667"
## CI_Gum     "0.000409035799509128" "Reject H0"         "0.0457990635262628"
## CD_Gum     "0.690840364244355"    "Fail to Reject H0" "0.839748890890375" 
##            [,4]                [,5]                [,6]               
## Kupiec_Gum "Reject H0"         "0.82296561431858"  "Fail to Reject H0"
## CI_Gum     "Reject H0"         "0.296111582252392" "Fail to Reject H0"
## CD_Gum     "Fail to Reject H0" "0.807570061778734" "Fail to Reject H0"

# Frank Copula
Kupiec_F<-c(test1_99_F$uc.LRp,test1_99_F$uc.Decision,test1_95_F$uc.LRp,test1_95_F$uc.Decision,test1_90_F$uc.LRp,test1_90_F$uc.Decision)
CI_F<-c(test1_99_F$cc.LRp,test1_99_F$cc.Decision,test1_95_F$cc.LRp,test1_95_F$cc.Decision,test1_90_F$cc.LRp,test1_90_F$cc.Decision)
CD_F<-c(test2_99_F$LRp,test2_99_F$Decision,test2_95_F$LRp,test2_95_F$Decision,test2_90_F$LRp,test2_90_F$Decision)

rbind(Kupiec_F,CI_F,CD_F)

##          [,1]                   [,2]                [,3]               
## Kupiec_F "8.1483655644754e-05"  "Reject H0"         "0.180919963625945"
## CI_F     "0.000409035799509128" "Reject H0"         "0.34583333112236" 
## CD_F     "0.438793457525467"    "Fail to Reject H0" "0.940312602513602"
##          [,4]                [,5]                [,6]               
## Kupiec_F "Fail to Reject H0" "0.323310244931084" "Fail to Reject H0"
## CI_F     "Fail to Reject H0" "0.238419945772615" "Fail to Reject H0"
## CD_F     "Fail to Reject H0" "0.98418169389875"  "Fail to Reject H0"

# Note: the expected number of exceedances is as following:
## VaR 99%: 17
## VaR 95%: 86
## VaR 90%: 173

exceed<-c(24,90,156,24,88,158,13,76,153,36,106,176,36,99,161)
Kupiec<-c(rep(0,9),1,1,0,1,0,0)
Christ_1<-c(rep(0,9),1,1,0,1,0,0)
Christ_2<-rep(0,15)
summary_backtesting<-cbind(exceed,Kupiec,Christ_1,Christ_2)
dimnames(summary_backtesting)<-list(c("Gaussian (VaR 99%)","Gaussian (VaR 95%)","Gaussian (VaR 90%)","T-Student (VaR 99%)","T-Student (VaR 95%)","T-Student (VaR 90%)","Clayton (VaR 99%)","Clayton (VaR 95%)","Clayton (VaR 90%)","Gumbel (VaR 99%)","Gumbel (VaR 95%)","Gumbel (VaR 90%)","Frank (VaR 99%)","Frank (VaR 95%)","Frank (VaR 90%)"),c("exceedance","UC","IND","CC"))
summary_backtesting

##                     exceedance UC IND CC
## Gaussian (VaR 99%)          24  0   0  0
## Gaussian (VaR 95%)          90  0   0  0
## Gaussian (VaR 90%)         156  0   0  0
## T-Student (VaR 99%)         24  0   0  0
## T-Student (VaR 95%)         88  0   0  0
## T-Student (VaR 90%)        158  0   0  0
## Clayton (VaR 99%)           13  0   0  0
## Clayton (VaR 95%)           76  0   0  0
## Clayton (VaR 90%)          153  0   0  0
## Gumbel (VaR 99%)            36  1   1  0
## Gumbel (VaR 95%)           106  1   1  0
## Gumbel (VaR 90%)           176  0   0  0
## Frank (VaR 99%)             36  1   1  0
## Frank (VaR 95%)             99  0   0  0
## Frank (VaR 90%)            161  0   0  0